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ABSTRACT 

This paper contains a critical comparison of estimators minimizing Wahba’s 
loss function. Some new results are presented for the QUaternion ESTimator 
(QUEST) and Estimators of the Optimal Quaternion (ESOQ and ESOQ2) to 
avoid the computational burden of sequential rotations in these algorithms. None 
of these methods is as robust in principle as Davenport’s q method or the 
Singular Value Decomposition (SVD) method, which are significantly slower. 

Robustness is only an issue for measurements with widely differing accuracies, 
so the fastest estimators, the modified ESOQ and ESOQ2, are well suited to 
sensors that track multiple stars with comparable accuracies. More robust forms 
of ESOQ and ESOQ2 are developed that are intermediate in speed. 

INTRODUCTION 

In many spacecraft attitude systems, the attitude observations are naturally represented as unit vectors. 
Typical examples are the unit vectors giving the direction to the sun or a star and the unit vector in the 
direction of the Earth’s magnetic Field. This paper will consider algorithms for estimating spacecraft attitude 
from vector measurements taken at a single time, which are known as “single-frame” methods or “point” 
methods, in contrast to filtering methods that employ information about spacecraft dynamics. Almost all 
single-frame algorithms are based on a problem proposed in 1965 by Grace Wahba 1 . Wahba’s problem is to 
find the orthogonal matrix A with determinant +1 that minimizes the loss function 

( 1 ) 

where [b ( } is a set of unit vectors measured in a spacecraft’s body frame, {rj are the corresponding unit 
vectors in a reference frame, and { a { ) are non-negative weights. In this paper we choose the weights to be 
inverse variances, a x - cr~ 2 , in order to relate Wahba’s problem to Maximum Likelihood Estimation 2 . This 
choice differs from that of Wahba and many other authors, who assumed the weights normalized to unity. 

It is possible and has proven very convenient to write the loss function as 

L(A) = A 0 -tr(A£ T ) 

with 

A o = S, a i 

and 

Now it is clear that L(A) is minimized when the trace, tr(AB T ), is maximized. 
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( 2 ) 

( 3 ) 

(4) 


This has a close relation to the orthogonal Procrustes problem, which is to find the orthogonal matrix A 
that is closest to B in the sense of the Frobenius (or Euclidean, or Schur, or Hilbert-Schmidt) norm 3 

MIf *£,>/,; =lr(MM T ). (5) 

Now 

\\A - B\\; = \\A\\l + 1 B||p - 2tr (AB r ) = 3 + ||B||~ - 2tr (AB J ) , (6) 

so Wahba’s problem is equivalent to the orthogonal Procrustes problem with the proviso that the 
determinant of A must be +1. 

The purpose of this paper is to give an overview in a unified notation of algorithms for solving Wahba’s 
problem, to provide accuracy and speed comparisons, and to present two significant enhancements of 
existing methods. The popular QUaternion Estimator (QUEST) and Estimators of the Optimal Quaternion 
(ESOQ and ESOQ2) algorithms avoid singularities by employing a rotated reference system. A method 
introduced in this paper uses the diagonal elements of the B matrix to determine a desirable reference 
system, avoiding expensive sequential computations. Also, tests show that a first-order expansion in the 
loss function is adequate, avoiding the need for iterative refinement of the loss function, and motivating the 
introduction of new first-order versions of ESOQ and ESOQ2, which are at present the fastest known first- 
order methods for solving Wahba’s problem. 

FIRST SOLUTIONS OF WAHBA’S PROBLEM 

J. L. Farrell and J. C. Stuelpnagel 4 , R. H. Wessner 5 , J. R. Velman 6 , J. E. Brock 7 , R. Desjardins, and 
Wahba presented the first solutions of Wahba’s problem. Farrell and Stuelpnagel noted that any real square 
matrix, including B, has the polar decomposition 

B = WR , (7) 

where W is orthogonal and R is symmetric and positive semidefinite. Then R can be diagonalized by 

R = VDV T y (8) 

where V is orthogonal and D is diagonal with elements arranged in decreasing order. The optimal attitude 
estimate is then given by 

A opt =Wdiag[l 1 det W]V T . (9) 

In most cases, detW is positive and A opt = W 7 but this is not guaranteed. Wessner proposed the solution 

A opt = (B T y\B T B) m = B(B T BT ia , (10) 

This has the disadvantage of requiring B to be non-singular, which means that a minimum of three vectors 
must be observed, although it is well known that two vectors are sufficient to determine the attitude. 

SINGULAR VALUE DECOMPOSITION (SVD) METHOD 

This method has not been widely used in practice, because of its computational expense, but it yields 
valuable analytic insights 8,9 . The matrix B has the Singular Value Decomposition 3 : 

B = U1V T = U diag[X n X 22 X 33 ] V T , (11) 

where U and V are orthogonal, and the singular values obey the inequalities X n > - £33 - 0. Then 

tr(/tfl T ) = tr(4 VdiagtS,, S 22 Z 33 ] f/ T ) = tr(t/ T A V'diagtl,, Z 22 S 33 ]). (12) 

The trace is maximized, consistent with the constraint det A = 1, by 

t/ r ^ op . ^ = diagfl 1 (det f/)(det V)] . (13) 

which gives the optimal attitude matrix 

A a p, = U diag[l 1 (det U)( det V)] V T . 


(14) 


The SVD solution is completely equivalent to the original solution by Farrell and Stuelpnagel, since 
Eq. (14) is identical to Eq. (9) with U = WV. The difference is that robust SVD algorithms exist now 3,10 . 

In fact, computing the SVD is one of the most robust numerical algorithms. 

It is convenient to define 

= X lh 5 2 = S 22 , and s 3 = (detLO(detV) X 33 , (15) 


so that > s 2 > |j 3 1 . We will loosely refer to s { , s 2 , and j 3 as the singular values, although the third singular 
value of B is actually |j 3 |. It is clear from Eq. (11) that redefinition of the basis vectors in the reference or 
body frame affects V or f/, respectively, but does not affect the singular values. 

The covariance of the rotation angle error vector in the body frame is given by 

P = U diag[(s 2 + j 3 )"’ (*3+5,)-’ (5,+j 2 r , lt/ T . (16) 

DAVENPORT’S q METHOD 

Davenport provided the first useful solution of Wahba’s problem for spacecraft attitude determination 11,12 . 

He parameterized the attitude matrix by a unit quaternion 13,14 


q l- 

esin( 0 / 2 ) 

q 4. 

cos (0 / 2 ) 


(17) 


as 


Mg) - (<?4 “ |q| 2 )7 + 2qq T — 2g 4 [q x] . 


(18) 


The rotation axis e and angle 0 will be useful later. Since A{q) is a homogenous quadratic function of < 7 , we 
can write 

tr(AS T ) = q T Kq (19) 

where K is the symmetric traceless matrix 


with 


K = 


S-ItrB 


z 

t vB 


S = B + B 


and 


z = 


5 23 B 32 

2.J 


= X i a i b . xr < 


( 20 ) 

( 21 ) 


It is easy to prove that the optimal unit quaternion is the normalized eigenvector of K with the largest 
eigenvalue, i.e, 9 the solution of 

* 9 .* •**.«*■ ( 22 ) 

With Eqs. (2) and (19), this gives the optimized loss function as 

UA 0 p ,) = A 0 -A mM . (23) 


Very robust algorithms exist to solve the symmetric eigenvalue problem 3,10 . 

The eigenvalues of the K matrix, A^ = A, > A 2 > A 3 > A 4 s A niiIl , are related to the singular values by 11 

Aj — S j "4" Sj 4" J' 3 , A 2 S j — S 2 — ^ 3 , A 3 — — 5j 4" ^2 J 3 , A 4 — — S | ^2 "4" 5' 3 . (24) 

The eigenvalues sum to zero because K is traceless. There is no unique solution if the two largest 
eigenvalues of K are equal, or + j 3 = 0. This is not a failure of the q method; it means that the data aren’t 
sufficient to determine the attitude uniquely. Equation (16) shows that the covariance is infinite in this case. 
This is expected, since the covariance should be infinite when the attitude is unobservable. 


QUATERNION ESTIMATOR (QUEST) 

This algorithm, first applied in the MAGSAT mission in 1979, has been the most widely used algorithm 
for Wahba’s problem 15,16 . Equation (22) is equivalent to the two equations 


and 


[(A lim +trB)/-S]q = < 7 4 z 
('Lux -trfl)<7 4 =q T z. 


Equation (25) gives 

q = + mi -S)- ] z = q 4 { adj [(A_ Hh tr B)I - S ]z}/det[( A max + tr 5)/ - 5] . 

The optimal quaternion is then given by 

*7opt 


1 

X 

<Jr 2+ H 2 

Y_ 


where 

and 

with 


X = adj [(A^ + tr B)I - 5] z = [al + ~ trfl) 5 + S 2 ] z 

Y s dettCA^ + tiB)I - 5] = + trfl) - det S , 

« = AL-(trB) 2 +tr(adjS). 


(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 


The second form on the right sides of Eqs. (29) and (30) follows from the Cayley-Hamilton Theorem 316 . 
These computations require knowledge of A^, which is obtained by substituting Eqs. (28) and (29) into 


Eq. (26), yielding: 


o = - trfl) - z T [a / + (A max - tr B) S + S 2 ] z . 


(32) 


Substituting a and y from Eqs. (30) and (31) gives a fourth-order equation in A^, which is simply the 
characteristic equation det(A max /- K) = 0. Shuster observed that A^ can be easily obtained by Newton- 
Raphson iteration of Eq. (32) starting from Aq as the initial estimate, since Eq. (23) shows that A max is very 
close to Aq if the optimized loss function is small 15,16 . In fact, a single iteration is generally sufficient. But 
numerical analysts know that solving the characteristic equation is one of the worst ways to find 
eigenvalues, in general, so QUEST is in principle less robust than Davenport’s original q method. The 
analytic solution of the characteristic equation is slower and no more accurate than the iterative solution. 

Shuster provided an estimate of the covariance of the rotation angle error vector in the body frame, 

(33) 

He also showed that 2 L(A opi ) obeys a j? probability distribution with 2n obi - 3 degrees of freedom, to a 
good approximation and assuming Gaussian measurement errors, where n obs is the number of vector 
observations 17 . This can often provide a useful data quality check. 

The optimal quaternion is not defined by Eq. (28) if y 2 + |x| 2 =0. Applying the Cayley-Hamilton theorem 
twice to eliminate S 4 and S 3 after substituting Eq. (29) gives, with some tedious algebra, 

y 2 +|x| 2 =y(dvr/dA), (34) 

where y/(A) is the quartic function defined implicitly by Eq. (32). It follows from the discussion following 
Eq. (15) that dyr / dA is invariant under rotations, since the coefficients in the polynomial V^) depend 
only on the singular values of 5, and dys I dA must be nonzero for the Newton-Raphson iteration for A^ 
to be successful 20 . The singular condition y 2 + |x| 2 = 0 is thus seen to be equivalent to y- 0, which means 
that (tf opt ) 4 = 0 and the optimal attitude represents a 180° rotation. Shuster devised the method of sequential 
rotations to avoid this singular case 15-18 . 


REFERENCE FRAME ROTATIONS 


The (<7 ) 4 = 0 singularity occurs because QUEST does not treat the four components of the quaternion on 
an equal basis. Davenport’s q method treats the four components symmetrically and is nonsingular, but 
some other methods have singularities similar to that in QUEST. These singularities can be avoided by 
solving for the attitude with respect to a reference coordinate frame rotated from the original frame by 180° 
about the x y y, or z coordinate axis. That is, we solve for one of the quaternions 


q ( ^q( 


0 


q 

<74 


e , 

0 


^e.-qxe, 

-q e, 


for i - 1, 2, 3, 


(35) 


where e, is the unit vector along the 7 th coordinate axis. We use the convention of Reference 14 for 
quaternion rotations, rather than the historic convention. These products are trivial to implement by merely 
permuting and changing signs of the of the quaternion components. For example, 

<7' =[<7i.<72><73’<74] T ®[1,0,0,0] t = [< 7 4 , - <7 3 , <7 2 , - <7,] t ■ (36) 


The equations for the inverse transformations are the same, since a 180° rotation in the opposite direction 
has the same effect. These rotations are also easy to implement on the input data, since a rotation about 
axis i y for example, simply changes the signs of the yth and kth columns of the B matrix, where {i,y, k } is 
a permutation of {1, 2, 3). 


The original QUEST implementation performed sequential rotations, one axis at a time, until an acceptable 
reference coordinate system was found. It is clearly preferable to avoid the computations required by 
sequential rotations by choosing a single desirable rotation as early in the computation as possible. The 
optimal rotated coordinate frame could be found by investigating the components of an a priori quaternion, 
which is always available in a star tracker application, since a fairly good a priori attitude estimate is needed 
to identify the stars in the tracker’s field of view. If the fourth component has the largest magnitude, no 
rotation is performed, while a rotation about the 7th axis is performed if the 7th component has the largest 
magnitude. Then Eq. (36) or its equivalent shows that the fourth component of the rotated quaternion will 
have the largest magnitude. This magnitude must be at least 1/2, but may not be any greater, because it is 
possible for all the components of a unit quaternion to have magnitude 1/2. 


The need for a priori attitude information can be avoided by substituting information from the B matrix. We 
can expand K in terms of its eigenvectors and eigenvalues as 


k = S Kwl = i £ *,< 4 *,*J 


i 












(37) 


where q M is the normalized eigenvector of K corresponding to the eigenvalue A^ and 


w 




^(<7^)23 ^(*7^)32 

^(<7^)31 ~ ^(< 7 /i )l 3 


(38) 


The second step in Eq. (37) follows from the fact that the eigenvalues sum to zero, and the third step 
follows from Eq. (18). We use Greek indices to label different quaternions, to avoid confusion with Latin 
indices used to label quaternion components. Comparison of Eqs. (37) and (20) establishes the identity 19 

8 = (39) 

From Eqs. (18) and (39), we have 


( 40 ) 


and 


25« =I>„M +trS ' ( 41 ) 

^ = 1 

These relations are interesting because the discussion following Eq. (15) makes it clear that the eigenvalues 
are not affected by a redefinition of the basis vectors in the body or reference frames, although the quaternion 
is affected by such a redefinition. Suppose now that we find the maximum of {5 U , fl 22 , Z? 33 ,trZ?} . If tr B is 
the maximum, and if is well separated from the other eigenvalues of the K matrix, Eqs. (40) and (41) 
suggest that the fourth component of q 0?t has the largest magnitude, so no rotation is necessary. If B lt is the 
maximum, the same considerations suggest that the ith component of q opi has the largest magnitude, so a 
rotation about the ith axis is performed. This will tend to put the largest component of q opt in the fourth 
position in the rotated frame, which means that the rotation angle in the rotated frame is small, and the 
180° rotation singularity is avoided. The reference system rotation is easily “undone” by Eq. (36) or its 
equivalent after the optimal quaternion has been computed. 

It will be useful for future reference to note that Eqs. (18), (21), and (39) give 

z= £■ ( 42 > 

M=1 


Using the orthonormality of the eigenvectors and Eq. (40), we find that 


and thus that 


M 2 =5X<^)4-(trfl) 2 . 

M=l 

H - mx|Aj * max(A imil -A min ). 


(43) 

(44) 


FAST OPTIMAL ATTITUDE MATRIX (FOAM) 

The SVD decomposition of B gives a convenient representation for adj£, detB, and ||£||p. These can be used 
to write the optimal attitude matrix as 20,21 

4* = (a*™ - det B)~ ] l(K + M 2 F )B + A^adj B J - BB J B] , (45) 

where 

K-iaL-Mp)- (46) 


It’s important to noie that all the quantities in Eqs. (45) and (46) can be computed without performing the 
SVD of B. In this method, A max is found from 

= tr (A op ,B T )=(KA^ - det + |K)||fif F + 3^^ det B - tr (Sfl T flfl T )], (47) 

or, after some matrix algebra, 

0 = </a mM ) = (AL -Ml) 1 -8A mM detB-4|adjS|| 2 F . (48) 

Equations (32) and (48) for V/(A inax ) would be numerically identical with infinite-precision computations, 
but the FOAM form of the coefficients is less subject to errors arising in finite-precision computations. 

The FOAM algorithm gives the convenient form for the error covariance: 

P = {Kk m - det BY\k1 + BB J ) . (49) 

A quaternion can be extracted from A opty with a cost of 13 MATLAB flops. This has two advantages: the 
four-component quaternion is often preferable to the nine-component attitude matrix, and the quaternion can 
be easily normalized if A opt is not exactly orthogonal due to computational errors 22 . 


ESTIMATOR OF THE OPTIMAL QUATERNION (ESOQ or ESOQ1) 

Davenport’s eigenvalue equation, Eq. (22), says that the optimal quaternion is orthogonal to all the 
columns of the matrix 

H = (50) 

which means that it must be orthogonal to the three-dimensional subspace spanned by the columns of H. 
The optimal quaternion is conveniently computed as the generalized four-dimensional cross-product of any 
three columns of this matrix 23 ' 25 . 

Another way of seeing this result is to examine the classical adjoint of H. Representing K in terms of its 
eigenvalues and eigenvectors and using the orthonormality of the eigenvectors gives, for any scalar A, 


adj( AT — A/) = adj 



= ^(A v -A)(A p -A)(A r -A)^^. 

/i=I 


( 51 ) 


where [pi, v , p, r] is a permutation of { 1, 2, 3, 4). Setting A = A miw = A, causes all the terms in this sum to 
vanish except the first, with the result 

adj H = a 2 ~ A_)(A 3 - A mM )(A 4 - (52) 

Thus g opt can be computed by normalizing any non-zero column of adj H , which we denote by index k. Let 
F denote the symmetric 3x3 matrix obtained by deleting the kih row and kih column from H, and let f 
denote the three-component column vector obtained by deleting the £th element from the kih column of H. 
Then the kih element of the optimal quaternion is given by 

07o P t )*=-cdetF, (53) 

and the other three elements are 

(tfopt)i,- -.jt-i.i+i, -,4 = c ( a dj/ r )f » (54) 

where the coefficient c is determined by normalizing the quaternion. It is desirable to let k denote the 
column with the maximum Euclidean norm, which is the column containing the maximum diagonal 
element of the adjoint, owing to the symmetry of H. Computing all the diagonal elements of adj//, though 
not as burdensome as QUEST’S sequential rotations, is somewhat expensive. It can be avoided by using the 
trace of the B matrix as in QUEST. In the ESOQ case, however, no rotation is performed; we merely 
choose k to be the index of the maximum element of {Z? n ,/? 22 , Z? 33 ,tr£} . 

The matrix F depends on the maximum eigenvalue A^; but it is interesting to note that f does not depend 
on A max , which only appears in the diagonal elements of H . The original formulation of ESOQ used the 
analytic solution of the characteristic equation 26 ; but the analytic formula sometimes gives complex 
eigenvalues, which is theoretically impossible for a real symmetric matrix. These errors arise from 
inaccurate values of the coefficients of the quartic characteristic equation, not from the solution method. It is 
faster, and equally accurate, to compute A^ by iterative solution of Eq. (48). Equation (32) would give a 
faster solution, but it would be less robust, and an even more efficient solution is described below. 

First Order Update (ESOQ1.1) 

Test results show that higher-order updates do not improve the performance of the iterative methods, 
providing motivation for developing a first-order approximation. The matrix H can be expanded to first order 
in AA s A 0 - as 

H = H° +(AA)/, (55) 

where 

H° = AT-A 0 /. (56) 

F= F° + (AA)/, 


Then 


(57) 


( 58 ) 


where F° is derived from H° in the same way that F is derived from H. Matrix identities give 

adj F = adj F° + AA[(tr F°)l-F°] 
and 

detF = detF 0 + (AA)tr(adjF 0 ), (59) 

to first order in A A . The characteristic equation can be expressed to the same order as 

0 = det H = (H° tt + AA) det F - f T (adj F) f = H° u detF 0 - f T g + [tf° tr(adjF°) + det F° - f T h] AA , (60) 

where 

g=(adjF°)f and h = [(trF°)I - F°]t. (61) 

Equation (60) is easily solved for AA = A 0 - A mM , and then the first order quaternion estimate is given by 

(<7o P ,)* = -c[detF° + (AA)tr(adjF 0 )] (62) 

and 

(^o P i ),,-,*-i,*+i.- .4 = c (g + h AA ) . (63) 

SECOND ESTIMATOR OF THE OPTIMAL QUATERNION (ESOQ2) 

This algorithm uses the rotation axis/angle form of the optimal quaternion, as given in Eq. (17). 
Substituting these into Eqs. (25) and (26) gives 

( A max - tr B) cos {<t> / 2) = z T e sin(0 / 2) (64) 

and 

z cos(0 / 2) = [( A^ + 1 tB)I - S] e sin (0 / 2) (65) 

Multiplying Eq. (65) by (A^ - tr B) and substituting Eq. (64) gives 

Mesin(0/2) = 0, (66) 

where 

M = (A maj - trB)[(A mM + tiB)I - S] - zz T = [m, • m 2 ■ m 3 ] . (67) 

These computations lose numerical significance if (A max - trB) and z are close to zero, which would be the 
case for zero rotation angle. We can always avoid this singular condition by using one of the sequential 
reference system rotations 15 " 18 to ensure that trB is less than or equal to zero. If we rotate the reference frame 
about the /th axis, then 

(tr 6) - B u -B m -B u = 2 B u - tr B, (68) 

where {/,y, k) is a permutation of { 1, 2, 3} and the components not marked “rotated” are defined with 
respect to the unrotated reference frame. Thus no rotation is performed if XrB is the minimum of 
{5 n ,5 22 ,B 3 3,trB} , while a rotation about the /th axis is performed if B a is the minimum. This will ensure 
the most negative value for the trace in the rotated frame. The rotation is easily “undone” by Eq. (36) or its 
equivalent after the quaternion has been computed. It follows from the orthogonality of the eigenvectors that 

Amin - trB ^ A^, (69) 

Equation (40) shows that trZ? = A^ if and only if q opt is the quaternion of the identity attitude, with 
cos(0/2) = 1, while tri? - A^ means that the identity quaternion is an eigenvector of K with eigenvalue 
Ay^jj,. The latter equality requires q opi to be the quaternion of a 180° rotation, with cos(0/2) = 0, since the 
orthogonality of the eigenvectors means that any two of the rotation matrices A(q M ) differ by a 180° 
rotation. However, q opt being the quaternion of a 180° rotation does not imply trB = A^ , but only the 
weaker condition A^ < tr B < A 2 . Thus a small value of trB is a better indicator of an attitude far from the 
identity than a large value of trB is of an attitude far from a 180° rotation, so this procedure should be even 
more reliable for avoiding singularities in ESOQ2 than in QUEST or ESOQ. 


Equation (66) says that the rotation axis is a null vector of M. The columns of adj M are the cross products 
of the columns of M: 

adj M = [m 2 xm 3 :m 3 xm t :m, x m 2 ]. (70) 

Because M is singular, all these columns are parallel, and all are parallel to the rotation axis e. Thus we set 

e = y/bi ( 71 ) 

where y is the column of adjM the cross product) with maximum norm. Because M is symmetric, it 
is only necessary to find the maximum diagonal element of its adjoint to determine which column to use. 
The rotation angle is found from Eq. (64) or one of the components of Eq. (65). Equations (40) and (44) 
show that choosing the rotated reference system that provides the most negative value of tri? makes Eq. (64) 
the best choice. With Eq. (71), this can be written 

( - tr B)|y| cos(0 / 2) = (z ■ y) sin(0 / 2) , (72) 

which means that there is some scalar rj for which 

cos (<j> / 2) = rj (z ■ y) (73) 

and 

sin(0 / 2) = 77 ( A^ - tr 5) |y| . (74) 

Substituting into Eq. (17) and using Eq. (71) gives the optimal quaternion as 27,28 



Note that it is not necessary to normalize the rotation axis. ESOQ2 does not define the rotation axis 
uniquely if M has rank less than two. This includes the usual case of unobservable attitude and also the case 
of zero rotation angle. Requiring tr B to be non-positive avoids zero rotation angle singularity, however. We 
compute by iterative solution of Eq. (48) in the general case, as for ESOQ. 

First Order Update (ESOQ2.1) 

The motivation for and development of this algorithm are similar to the those of the first order update used 
in ESOQ 1.1. The matrix M can be expanded to first order in AA = A 0 - A max as 

M = M° + (AA)AT, (76) 

where 

M° = (A 0 - tr5)[(A 0 4 - tiB)I - S] - zz T = [m? : m° 2 : mj] (77) 

and 

N = S-2A 0 / = [n i :n 2 :n 3 ]. (78) 

To this same order, we have 

y s m t x = (m° + n ( AA) x (mj + n J AA) = y° + p AA , (79) 

where 

y°=mj ) xinj and psmfxn.+n^xmj. (80) 

The maximum eigenvalue can be found from condition that M is singular; to first-order: 

0 = det M = (m J xm y ) m k = (y° +pAA)-(m° +n*AA) = y° -m° + (y° n k +m° k *p)AA, (81) 

where {/,y, k} is a cyclic permutation of { 1, 2, 3}. This is solved for AA s A 0 - A^, and then the attitude 
estimate is found by substituting Eq. (79) and A^ = A 0 - AA into Eq. (75). 



There is an interesting relation between the eigenvalue condition det A/(A) = 0 used in ESOQ2. 1 and the 
condition y/(A) = 0 used in other algorithms. Straightforward matrix algebra shows 

detM(A) = (A-tr£)V(A). ( 82 ) 


Thus det /V/(A) has the four roots of i/r(A), the eigenvalues of Davenport’s K matrix, and an additional 
double root at trS. Equation (40) shows that trZ? depends on the reference frame axes, and choosing the 
reference axes maximizing -tr B assures that these two spurious roots are far from the desired root at A^. 

TESTS 

We test the accuracy and speed of MATLAB implementations of these methods, using simulated data. The q 
and SVD methods use the functions eig and svd, respectively; the others use the equations in this paper. 
MATLAB uses IEEE double-precision floating-point arithmetic, in which the numbers have approximately 
16 significant decimal digits 29 . 

We analyze three test scenarios. In all these scenarios, the pointing of one spacecraft axis, which we take to 
be the spacecraft x axis, is much better determined that the rotation about this axis. This is a very common 
case that arises in spacecraft that point a single instrument (like an astronomical telescope) very precisely. 
This is also a characteristic of attitude estimates from a single narrow-field-of-view star tracker, where the 
rotation about the tracker boresight is much less well determined than the pointing of the boresight. The x 
axis error and the yz error, which is the error about an axis orthogonal to the x axis and determines the x 
axis pointing, are computed from an error quaternion q cn by writing 


'q.rr 1 = ["«* sin(0, / 2)1 r sin(0 yz / 2) 
q, ni \ L cos(0, 12) \ [ cos(0 >; / 2) _ 

"e, cos (0 >z / 2)sin(0 y / 2) + e yz cos (0, / 2 ) 8111 ( 0 ^, / 2) - (e, x e yz )sin(0 x / 2)sin(0 y _ / 2) 

cos (<j) x /2)cos(0^ 12) 


where e, = [1 0 0] T and e yz is a unit vector orthogonal to e,. We can always find <p x in [- n , n] and </> yz in 
[0, 7i ] by selecting e yz appropriately. Then, since and x e yz form an orthonormal basis for the y-z 
(or 2-3) plane, the error angles are given by 

<t> I =2Un-'(q crrl /q ari ) (84) 


and 



(85) 


Equations (84) and (85) would be unchanged if the order of the rotations about e x and e yz were reversed; only 
the unit vector would be different. 

The total error, which is the principal angle of the rotation represented by the error quaternion, is given by 

COS (0^ / 2) = q^ t = cos (0, / 2) cos(0 y , / 2) . (86) 

Thus 0^, / 2 is the hypotenuse of a right spherical triangle with sides <j> x / 2 and 0^/2 .This is the 
spherical trigonometry equivalent of taking two orthogonal components of an error vector. 

First Test Scenario 

The first scenario simulates a single star tracker with a narrow field of view and boresight at [1, 0, 0] r . This 
is an application for which the QUEST algorithm has been widely used. We assume that the tracker is 
tracking five stars at 

'll [0.99712] [ 0.99712] [0.99712] [ 0.99712" 

b, = 0 , b 2 = 0.07584 , b 3 = -0.07584 , b 4 = 0 , and b 5 = 0 . (87) 

oj |_ o J [ o J |_°- 07584 J L-°- 07584 _ 


We simulate 1000 test cases with uniformly distributed random attitude matrices, which we use to map the 
five observation vectors to the reference frame. We add Gaussian random noise with equal standard 
deviations of 6 arcseconds per axis to the reference vectors rather than the observation vectors, so that Eq. 
(87) will remain valid in the presence of noise, and then normalize the reference vectors. 

The loss function is computed with measurement variances in (radians) 2 , since this results in 2L{A 0?K ) 
approximately obeying a distribution. The minimum and maximum values of the loss function in the 
1000 test runs are 0.23 and 12.1, respectively. The probability distribution of the loss function is plotted as 
the solid line in Figure 1, and several values of P(x 2 | v) for X 2 = 2L(A opt ) and v = 2n obs -3 = 7 are 

plotted as circles 26 . The agreement is seen to be excellent. 

The estimation errors in arcseconds for the star tracker scenario are presented in Table 1, as both the RSS 
(outside of parentheses) and the maximum (in parentheses) over the 1000 cases. The q method and the SVD 
method should both give the truly optimal solution, since they are based on robust matrix analysis 
algorithms 3,10 . The q method is taken as optimal by definition, so no estimated-to-optimal differences are 
presented for that algorithm, and the differences between the SVD and q methods provide an estimate of the 
computational errors of both methods. In particular, the loss function is computed exactly by both methods, 
in principle, which means in practice that it is computed to about one part in 10 5 . Equation (3) gives 
A 0 = 5.9 x 10 9 rad' 2 for this scenario, so this is the expected accuracy in double-precision machine 
computations. No estimate of the loss function is provided when no update of is performed, 
accounting for the lack of entries in the loss function column in the tables for these cases. 



Figure 1: Empirical (solid line) and Theoretical (dots) Loss Function Distribution 
for the Seven-Degree-of-Freedom Star Tracker Scenario 


Not all the decimal places exhibited are significant, since the results of 1000 different random cases would 
not agree with these to more than two decimal places. The extra decimal places are shown to emphasize the 
fact that although the different algorithms give results that are closer or farther from the optimal estimate, 
all the algorithms provide estimates that are equally close to the true attitude. The differences between the 
estimated and optimal values further show that one Newton-Raphson iteration for A, nax is always sufficient; 
a second iteration provides no practical improvement in the estimate for this scenario. 

Equation (33) gives the covariance for the star tracker scenario as 

5 

P = (6 arcsec) 2 [5/ - J^b.-bT]” 1 = diag[1565,7.2,7.2]arcsec 2 , (88) 

/=i 

which gives the standard deviations of the attitude estimation errors as 

a x = V 1565 arcsec = 40 arcsec and o yz = V7.2 + 7.2 arcsec = 3.8 arcsec . (89) 


It is apparent that this covariance estimate is quite accurate. 

Second Test Scenario 


The second scenario uses three observations with widely varying accuracies to provide a difficult test case 
for the algorithms under consideration. The three observation vectors are 


T 


'-0.99712' 


'-0.99712' 

0 

, b 2 

0.07584 

, and b 3 = 

-0.07584 

_o_ 


0 


0 


(90) 


Table 1: Estimation Errors for Star Tracker Scenario 


Algorithm 

iterations 

RSS (max) estimated-to-optimal 

RSS (max) estimated-to-true 

loss function 

x (arcsec) 

yz (arcsec) 

x (arcsec) 

yz (arcsec) 

Q 

— 

— 

— 

— 

38.41 (122.5) 

3.829 (9.252) 

SVD 

— 

0.4 (1.8) x 10~ 5 

1.4 (5.6) x lO" 8 

0.8 (2.9) x lO" 10 

38.41 (122.5) 

3.829 (9.252) 

FOAM 

0 

— 

0.014 (0.078) 

9.7 (36) x lO’ 5 

38.41 (122.5) 

3.829 (9.252) 

1 

0.4 (1.6) x lO' 5 

1.5 (5.6) x lO' 8 

26 (104) x 10-‘° 

38.41 (122.5) 

3.829 (9.252) 

2 

0.4 (1.7) x lO' 5 

1.5 (5.6) x lO' 8 

26 (88) x 10-'° 

38.41 (122.5) 

3.829 (9.252) 

QUEST 

0 

— 

0.015 (0.078) 

9.6 (36) x lO" 5 

38.41 (122.5) 

3.829 (9.252) 

1 

2.5 (7.2) x 10~ 5 

10.1 (46) x 10~ 8 

6.1 (26) x lO- 10 

38.41 (122.5) 

3.829 (9.252) 

2 

2.9 (8.4) x 10~ 5 

11.1 (50) x 10' 8 

7.2 (25) x 10-'° 

38.41 (122.5) 

3.829 (9.252) 

ESOQ2 

0 

— 

0.014 (0.078) 

9.7 (36) x 10 5 

38.41 (122.5) 

3.829 (9.252) 

1 

0.4 (1.7) x 10~ 5 

1.5 (6.1) x lO' 8 

2.0 (10) x 10-'° 

38.41 (122.5) 

3.829 (9.252) 

2 

0.4 (1.8) x 10‘ 5 

1.5 (6.1) x lO’ 8 

2.1 (ll)x lO' 10 

38.41 (122.5) 

3.829 (9.252) 

ESOQ2.1 

1 

0.4 (1.6) x lO' 5 

1.5 (5.9) x lO' 8 

1.9 (12) x 10-‘° 

38.41 (122.5) 

3.829 (9.252) 

ESOQ 

0 

— 

0.015 (0.078) 

9.6 (36) x 10* s 

38.41 (122.5) 

3.829 (9.252) 

1 

0.4 (1.7) x 10~ 5 

1.5 (6.2) x 10- 8 

9.6 (39) x lO’ 10 

38.41 (122.5) 

3.829 (9.252) 

2 

0.4 (1.8) x 10- 5 

1.5 (6.2) x lO’ 8 

9.6 (39) x lO 10 

38.41 (122.5) 

3.829 (9.252) 

ESOQ1.1 

1 

1.0 (5.3) x 10~ ! 

4.1 (24) x 10‘ s 

7.0 (29) x lO’ 10 

38.41 (122.5) 

3.829 (9.252) 









































































































We simulate 1000 test cases as in the star tracker scenario, but with Gaussian noise of one arcsecond per 
axis on the first observation, and 1° per axis on the other two. This models the case that the first 
observation is from an onboard astronomical telescope, and the other two observations are from a coarse sun 
sensor and a magnetometer, for example. A very accurate estimate of the orientation of the x axis is required 
in such an application, but the rotation about this axis expected to be fairly poorly determined. 

The minimum and maximum values of the loss function computed by the q method in the 1000 test runs 
for the second scenario are 0.003 and 8.5, respectively. The probability distribution of the loss function is 
plotted as the solid line in Figure 2, and several values of the x 2 distribution with three degrees of freedom 
are plotted as circles. The agreement is almost as good as the seven-degree-of- freedom case. 

The estimation errors for this scenario are presented in Table 2, which is similar to Table 1 except that the 
rotation errors about the x axis are given in degrees. The agreement of the q and SVD methods is virtually 
identical to their agreement in the star tracker scenario, but the other algorithms show varying performance. 
Equation (3) gives A 0 = 8.5 x 10 10 rad -2 for this scenario, so the expected accuracy of the loss function in 
double-precision machine computations is on the order of 10' 5 , which is the level of agreement between the 
values computed by the q and SVD methods. None of the other methods computes the loss function nearly 
as accurately. This differs from the first scenario, where all the algorithms came close to achieving the 
maximum precision available in double-precision arithmetic. 



Figure 2: Empirical (solid line) and Theoretical (dots) Loss Function Distribution 
for the Three-Degree-of-Freedom Unequal Measurement Weight Scenario 




The iterative computation of A max in QUEST, ESOQl.l, and ESOQ2.1 is poor, but this has surprisingly 
little effect on the determination of the x axis pointing. The determination of the rotation about the x axis 
is adversely affected by an inaccurate computation of A, ntu , however, with maximum deviations from the 
optimal estimate of almost 180°. The only useful results of QUEST are obtained by not performing any 
iterations for A mjx . The iterative computation of A max by Eq. (48) in FOAM, ESOQ, and ESOQ2 improves 
the agreement with the optimal estimate, but does not result in better agreement with the true attitude. It 
seems probable that Eq. (48) provides a better estimate of A, IUIX because it deals with B directly, while the 
other algorithms use the symmetric and skew parts S and z instead. 

The predicted covariance in this scenario is, to a very good approximation, 

P = diag[A(l - 0.997 12 2 )" 1 deg 2 ,l arcsec 2 ,! arcsec 2 ], (91) 

which gives 

a x = 9.3 deg and a yz =1.4 arcsec , (92) 

in agreement with the best results in Table 2. 

Third Test Scenario 

The third scenario investigates the effect of measurement noise mismodeling, illustrating problems that first 
appeared in analyzing data from the Upper Atmosphere Research Satellite 30 . Of course, no one would 
intentionally use erroneous models, but it can be very difficult to determine an accurate noise model for real 
data, and the assumption of any level of white noise is often a poor approximation to real measurement 
errors. This scenario uses the same three observation vectors as the second scenario, given by Eq. (90). We 
again simulate 1000 test cases, but with Gaussian white noise of 1° per axis on the first observation and 
0.1° per axis on the other two. The estimator, however, incorrectly assumes measurement errors of 0.1° per 
axis on all three observations, so it weights the measurements equally. 


Table 2: Estimation Errors for Unequal Measurement Weight Scenario 


Algorithm 
Ana* iterations 

RSS (max) estimated-to-optimal 

RSS (max) estimated-to-true 

loss function 

* (deg) 

yz (arcsec) 

X (deg) 

yz (arcsec) 

<7 

- 

— 

— 

— 

9.5 (34) 

1.42 (3.57) 

SVD 

- 

1.6 (6.9) x 10“ 5 

1.4 (8.0) x 10~ 5 

7.7 (24) x 10-“ 

9.5 (34) 

1.42 (3.57) 

FOAM 

0 

— 

1.5 (9.9) 

7.9 (29) x 10 3 

9.5 (34) 

1.42 (3.57) 

1 

0.09 (0.7) 

0.09 (1.1) 

8.0 (26) x lO' 3 

9.5 (34) 

1.42 (3.57) 

2 

0.0007 (0.012) 

0.0008 (0.013) 

7.8 (29) x lO' 3 

9.5 (34) 

1.42 (3.57) 

QUEST 

0 

— 

1.9(12) 

0.4 (3.8) x 1(T 3 

9.6 (34) 

1.42 (3.57) 

1 

768 (2329) 

60 (170) 

2.3 (9.0) x lO' 3 

48 (90) 

1.42 (3.57) 

2 

1796 (38501) 

62(175) 

4.8 (95) x 10' 3 

48 (91) 

1.42 (3.57) 

ESOQ2 

0 

— 

1.5 (9.9) 

1.1 (6.8) x lO' 3 

9.5 (34) 

1.42 (3.57) 

1 

0.09 (0.7) 

0.09 (1.1) 

1.3 (9.4) x lO' 3 

9.5 (34) 

1.42 (3.57) 

2 

0.0007 (0.012) 

0.0008 (0.013) 

1.1 (7.1) x lO' 3 

9.5 (34) 

1.42 (3.57) 

ESOQ2.1 

1 

59 (370) 

39(178) 

1.5 (14) x 10~ 3 

29 (91) 

1.42 (3.57) 

ESOQ 

0 

— 

1.9(12) 

4.8 (23) x lO' 3 

9.6 (34) 

1.42 (3.57) 

1 

0.09 (0.7) 

0.10(1.1) 

5.3 (28) x lO' 3 

9.5 (34) 

1.42 (3.57) 

2 

0.0007 (0.012) 

0.0008 (0.013) 

5.2 (24) x 10' 3 

9.5 (34) 

1.42 (3.57) 

ESOQl.l 

1 

327(1727) 

60(177) 

2.6 (34) x lO' 3 

43 (90) 

1.42 (3.57) 



































































































The minimum and maximum values of the loss function computed by the q method in the 1000 test runs 
for the third scenario are 0.07 and 453, respectively. The probability distribution of the loss function is 
plotted in Figure 3. The theoretical three-degree-of- freedom distribution is not plotted, since it would be a 
very poor fit to the data. More than 95% of the values of L(A opt ) are theoretically expected to lie below 4, 
according to the % l distribution plotted in Figure 2, but almost half of the values of the loss function 
plotted in Figure 3 have values greater than 50. Shuster has emphasized that large values of the loss 
function are an excellent indication of measurement mismodeling or simply of bad data. 

The estimation errors for this scenario are presented in Table 3, which is similar to Tables 1 and 2 except 
that all the angular errors are given in degrees. The truly optimal q and SVD methods agree even more 
closely than in the other scenarios. Equation (3) gives A 0 = 5 x 10 5 rad" 2 for this scenario, so the expected 
accuracy of the loss function in double-precision machine computations is on the order of 10" l °, the level of 
agreement between the q and SVD methods. As in the second scenario, none of the other methods computes 
the loss function nearly as accurately. In the third scenario, though, the iterative computation of X^ works 
well for all the algorithms, and both iterations improve the agreement of the loss function and attitude 
estimates with the optimal values. The first order refinement is reflected in a reduction of the attitude errors, 
particularly in determining the rotation about the x axis, but no algorithm is aided significantly by a 
second-order update. As in the first scenario, all the algorithms with the first order update to X^ perform as 
well as the q and SVD methods. 



Figure 3: Empirical Loss Function Distribution for the 
Mismodeled Measurement Weight Scenario 




Speed 

There are two caveats to make with regard to timing comparisons. First, absolute speed numbers are not 
very Important for ground computations, since the actual estimation algorithm is only a part of the overall 
attitude determination data processing effort. Absolute speed was more important in the past, when 
thousands of attitude solutions had to be computed by slower machines, which is why QUEST was so 
important for the MAGSAT mission. Second, the longest computation time is more important than the 
average time for a real-time computer in a spacecraft attitude control system or a star tracker, which must 
finish all its required tasks in a limited time. This penalizes QUEST for real-time applications, unless we 
use an a priori attitude estimate or information from the B matrix to eliminate the need for sequential 
rotations, as described above. 

Figures 4 and 5 show the maximum number of MATLAB floating-point operations (flops) to compute an 
attitude using two to six reference vectors; the times to process more than six vectors follow the trends seen 
in the figure. The inputs for the timing tests are the n obs normalized reference and observation vector pairs 
and their n obs weights. One thousand test cases with random attitudes and random reference vectors with 
Gaussian measurement noise were simulated for each number of reference vectors. 

Figure 4 plots the times of the more robust methods. The break in the plots for FOAM, ESOQ, and 
ESOQ2 at n obs = 3 results from using an exact solution of the characteristic equation in the two-observation 
case, when det B = 0 and Eq. (48) shows that l/r(A max ) is a quadratic function of For three or more 
observations, these algorithms are timed for a first-order update to A^ using Eq. (48). Additional iterations 
for A nux are not expensive, however, costing only 1 1 flops each. It is clear that the q method and the S VD 
method require significantly more computational effort than the other algorithms, as expected. The q method 
is more efficient than the SVD method, except in the least interesting two-observation case. The other three 
algorithms are much faster, with the fastest, ESOQ and ESOQ2, being nearly equal in speed. 


Table 3: Estimation Errors for Mismodeled Measurement Weight Scenario 


Algorithm 
A^ iterations 

RSS (max) estimated-to-optimal 

RSS (max) estimated-to-true 

loss function 

x (deg) 

yz (deg) 

x (deg) 

yz (deg) 

<7 

— 

— 

— 

— 

0.96 (3.62) 

0.49 (1.17) 

SVD 

— 

4.1 (22) x 1CT 10 

3.8 (17) x 10- 12 

2.3 (7.3) x 10- 14 

0.96 (3.62) 

0.49 (1.17) 

FOAM 

0 

— 

0.7 (5.9) 

4.0 (21) x 10" 3 

1.18 (5.42) 

0.49 (1.16) 

1 

2.6 (24) 

0.020 (0.33) 

1.0 (11) x 10" 4 

0.96 (3.60) 

0.49 (1.17) 

2 

0.004 (0.07) 

0.4 (10) x 10- 4 

1.7 (35) x 10" 7 

0.96 (3.62) 

0.49 (1.17) 

QUEST 

0 

— 

0.9 (7.6) 

4.3 (29) x 10‘ 3 

1.27 (7.66) 

0.49 (1.16) 

1 

2.6 (24) 

0.023 (0.33) 

1.1 (11) x 10* 4 

0.96 (3.60) 

0.49 (1.17) 

2 

0.004 (0.07) 

0.4 (10) x 10- 4 

1.7 (35) x lO' 7 

0.96 (3.62) 

0.49(1.17) 

ESOQ2 

0 

— 

0.7 (5.9) 

4.0 (21) x 10~ 3 

1.18 (5.42) 

0.49 (1.16) 

1 

2.6 (24) 

0.020 (0.33) 

1.0 (11) x 10" 4 

0.96 (3.60) 

0.49 (1.17) 

2 

0.004 (0.07) 

0.4 (10) x lO" 4 

1.7 (35) x lO' 7 

0.96 (3.62) 

0.49 (1.17) 

ESOQ2.1 

i [ 

2.6 (24) 

0.020 (0.33) 

0.6 (5.8) x 10- 4 

0.96 (3.60) 

0.49 (1.17) 

ESOQ 

0 

— 

0.9 (7.6) 

4.3 (29) x 10- 3 

1.27 (7.66) 

0.49 (1.16) 

1 

2.6 (24) 

0.023 (0.33) 

1.1 (11) x lO" 4 

0.96 (3.60) 

0.49 (1.17) 

2 

0.004 (0.07) 

0.4 (10) x 10" 4 

1.7 (35) x 10" 7 

0.96 (3.62) 

0.49 (1.17) 

ESOQ 1.1 

1 

2.6 (24) 

0.023 (0.33) 

1.3 (27) x 10-* 

0.96 (3.60) 

0.49 (1.17) 










































































































Figure 5 compares the timing of the fastest methods, which generally use zeroth order and First order 
approximations for A llllw . Both QUEST(i) and ESOQ2.1 use the exact quadratic solution for A inax in the 
two-observation case, but ESOQ 1.1 uses its faster first order approximation for any number of 
observations. It is clear that ESOQ and ESOQ2 are the fastest algorithms using the zeroth order 
approximation for A niiU , and ESOQ 1.1 is the fastest of the first order methods. 

CONCLUSIONS 

This paper has examined the most useful algorithms for estimating spacecraft attitude from vector 
measurements based on minimizing Wahba’s loss function. These were tested in three scenarios, which 
show that the most robust, reliable, and accurate estimators are Davenport’s q method and the Singular 
Value Decomposition (SVD) method. This is not surprising, since these methods are based on robust and 
well-tested general-purpose matrix algorithms. The q method, which computes the optimal quaternion as the 
eigenvector of a symmetric 4x4 matrix with the largest eigenvalue, is the faster of these two. 

Several algorithms are significantly less burdensome computationally than the q and SVD methods. These 
methods are less robust in principle, since they solve the quartic characteristic polynomial equation for the 
maximum eigenvalue, a procedure that is potentially numerically unreliable. Algorithms that use the form 
of the characteristic polynomial from the Fast Optimal Attitude Matrix (FOAM) algorithm performed as 
well as the q and SVD methods in practice, however. The fastest of these algorithms are the Estimators of 
the Optimal Quaternion, ESOQ and ESOQ2. 



Figure 4: Execution Times for Robust Estimation Algorithms 

FOAM, ESOQ, and ESOQ2 use first order approximation for 



Part of the speed advantage of the methods tested in this paper over previous studies results from using the 
information contained in the diagonal elements of B to eliminate sequential rotations in QUEST and extra 
computations in ESOQ and ESOQ2. It would be equally fast, and more accurate in principle, to use an a 
priori quaternion to find the best rotated reference frame. Simulations comparing the methods using the 
diagonal elements of B with methods using information contained in an a priori quaternion showed equally 
accurate estimates with the two approaches, however, and the former method is preferable for its generality. 

All the algorithms tested perform as well as the more robust algorithms in cases where measurement 
weights do not vary too widely and are reasonably well modeled. These include most of the cases for which 
vector observations are used to compute spacecraft attitude, in particular the case of an attitude solution 
from multiple stars. If the measurement uncertainties are not well represented by white noise, however, an 
update is required, while this update can be disastrous if the measurement weights span a wide range. 

The examples in the paper show that these robustness concerns are not an issue for the processing multiple 
star observations with comparable accuracies, the most common application of Wahba’s loss function. 

Thus the fastest algorithms, the zeroth-order ESOQ and ESOQ2 and the first-order ESOQ 1. 1, are well suited 
to star tracker attitude determination applications. In general-purpose applications where measurement 
weights may vary greatly, one of the more robust algorithms may be preferred. 



number of observed vectors 


Figure 5: Execution Times for Fast Estimation Algorithms 

Numbers in parentheses denote order of approximation. 
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